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In recent years, many chemical reactions have been studied by means of the quasi- 
Qh' classical trajectory (QCT) method within the Gaussian binning (GB) procedure. 

The latter consists in "quantizing" the final vibrational actions in Bohr spirit by 
putting strong emphasis on the trajectories reaching the products with vibrational 
actions close to integer values. A major drawback of this procedure is that if iV 
is the number of product vibrational modes, the amount of trajectories necessary 
to converge the calculations is ~ 10^ larger than with the standard QCT method. 



Applying it to polyatomic processes is thus problematic. In a recent paper, however, 



Czako and Bowman propose to quantize the total vibrational energy instead of the 
vibrational actions [G. Czako and J. M. Bowman, J. Chem. Phys., 131, 244302 
(2009)], a procedure called 1GB here. The calculations are then only ~ 10 times 
more time-consuming than with the standard QCT method, allowing thereby for 



considerable numerical saving. In this paper, we propose some theoretical arguments 
supporting the 1GB procedure and check its validity on model test cases as well as 



the prototype four-atom reaction OH+D2 — > HOD+D. 



I. INTRODUCTION 

Improving our ability to accurately describe gas-phase chemical reactions and inelastic 
collisions is a stimulating theoretical issue at the interface of physics and chemistry |l| and 
a necessary step towards a deep understanding of the evolution of planetary atmospheres 
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and interstellar clouds. 

Assuming that for a given process, the electronic problem has been solved |2|, i.e., the 
potential energy of interaction between nuclei is known, nuclear motions can be studied 



either quantum 



Q r 



3HlO| or classical mechanically [111 Il2| . For the present time, however, 
quantum scattering approaches can hardly be applied to more than three-atom processes, 
despite current computer performances and a great deal of methodological effort made to 



go beyond the triatomic problem 



13 
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On the other hand, the classical approach, well known as the quasi-classical trajectory 
method (QCTM) |lll.ll2j|. is much less time consuming and can therefore be applied to almost 
any process, independently on the number of atoms involved. We focus our attention on 
this method in the present paper. 

A major goal of QCTM is to predict the distributions of the translational energy between 
bimolecular collision or photodissociation products as well as the distribution of their quan- 
tum states [1]. These distributions, measured in molecular beam experiments, are among 
the most refined data on chemical reactivity and molecular reaction dynamics. In this work, 
we concentrate on the possible descriptions of these two quantities within QCTM. 

In its standard implementation, QCTM deals with the standard binning (SB) procedure 
(or histogram method) for assigning trajectories to the various quantum states available. In 
order to introduce this procedure, we consider the three-atom exchange reaction of the type 
A + BC — > AB + C. If at the end of a given reactive trajectory, the vibrational action of 
AB is x in units of h (see appendix A for the mathematical definition of x) and its rotational 
angular momentum is j in units of H, the trajectory is assumed to only contribute to the 
AB quantum state (x,j) where x and j are the nearest integers of x and j respectively (in 
the following, the nearest integer of any variable will also be denoted by the variable with a 
bar on top of it). 

About ten years ago, however, it was suggested that such a procedure might lead to 
wrong predictions when the energy available to the products is too low for the quantum and 
classical densities of product states to be equal, or equivalently, when the available quantum 
states are widely spaced as compared to the energy disposal [18]. A Gaussian Binning 
(GB) procedure was then proposed [l8j| which amounts to assigning to each trajectory a 
Gaussian statistical weight such that the closer the final actions to their nearest integers, 
the larger the weight (by action, we mean here both vibrational actions and rotational 
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angular momenta in the previously defined units). For the previous triatomic process, the 
Gaussian weight of the trajectory ending with (x,j) is 

G(x,j)=G(x-x)GV-j) (1) 

with 



G ( u ) = ^7' ( 2 ) 



ir^'e 



e being usually kept at ~ 0.05 |l9l-|2l|. Like in the SB procedure, trajectories do only 
contribute to the quantum state defined by the center (x,j) of the bin or unit square in 
which (x,j) stands. The GB procedure is therefore a practical way of taking into account 
Bohr quantization in the analysis of the final results. The GB procedure turns out to be 



22|, a 



a reminiscence of the use of narrow boxes proposed by Ron et al in the early 80' s 
method apparently ignored or forgotten by QCTM users. 

Initially proposed on the basis of intuitive arguments, the GB procedure was later shown 
to be a practical implement at ion of classical S matrix theory (CSMT) in the random phase 



approximation 
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24], CSMT being the first and simplest (or least complex) semi-classical 



app roach of molecular collisions pioneered by Miller and Marcus in the early seventies 



3l|. 
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The Gaussian weight G(u) is characterized by a full width at half maximum of ~ 10 
percent. This means that the values of x and j respectively in the ranges [x — 0.05, x + 0.05] 
and [j — 0.05, j + 0.05] mostly contribute to the GB population of the level (x, j), as compared 
with the values in the unit ranges [x — 0.5, x + 0.5] and [j — 0.5, j + 0.5] which contribute to 
the SB population. Therefore, the area in the (x, j) plane contributing to the GB population 
is ~ 100 times smaller than the one contributing to the SB population and it is necessary 
to run ~ 100 times more trajectories within the GB procedure than within the SB one for 
the same level of convergence of the final results. 

In many experiments, however, the number of available rotational states of AB is signif- 
icantly larger than the number of its vibrational states (more than ~ 10 against less than 
~ 3) and one arrives at the same result when weighting the trajectories by Eq. ([1]) or by 
G(x — x) alone. Within this partial GB procedure, corresponding to Eqs. (13) and (14) of 
reference 24, it is thus sufficient to run ~ 10 times more trajectories than within the HB 
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one |l9|-|2l|,y,|32|-|35|. 



However, considering polyatomic reactions where the number of vibrational modes is 
easily ten or more, strongly clouds the situation. The reason is that "quantizing" N modes 
amounts to weight the trajectories by a product of iV Gaussians. Therefore, one is led to 
run ~ 10^ times more trajectories within the GB procedure than within the SB one. For 
the reaction F+CFL — > FH+CH3 and its isotopic variants, much studied experimentally 
in the recent years [36], the previous number is ~ 10 millions ! Since one needs at least a 
few hundreds of thousands of trajectories within the HB procedure, one should run a few 
trillions of trajectories within the GB procedure, which is just not feasible. 

In order to circumvent this difficulty, Czako and Bowman recently proposed to weight 
the trajectories by G{u) (see Eq. ([2])) with 

X{ being the vibrational action for the i th mode and a;* the corresponding frequency [37|. In 
other words, they proposed to quantize, with one Gaussian only, the total vibrational energy 
(in the harmonic approximation) instead of the vibrational actions. Consequently, this 1GB 
procedure allows for a huge amount of computational savings for large systems. 

The goal of the present paper is to propose theoretical arguments supporting this proce- 
dure and check its validity on model as well as actual processes. 

The paper is organized as follows. In section II, the 1GB procedure is shown to be 
equivalent to the usual GB procedure for statistical collinear processes. We then discuss 
the conditions for its validity in the general case. The predictions to which it leads are 
compared in section III with the usual SB and GB predictions for a model test case involving 
three vibrational modes. In section IV, the approach is applied to the prototype four-body 
chemical reaction OH+D 2 — > HOD+D which is among the simplest polyatomic bimolecular 



reactions 



38 



44] . We finally conclude in section V. 



II. THEORETICAL ANALYSIS OF THE 1GB PROCEDURE 

In a first step, we focus our attention on collinear processes in the course of which nuclei 
keep on a line fixed in the laboratory frame. The realistic three-dimensional case where 



A Collisional system invoHin^tiEOMEmeiM lAdMLYSIS OF THE 1GB PROCEDURE 
rotation motions are active is considered in a second step. 

A. Collisional system involving two vibrational modes 

Consider the collinear inelastic collision between atom A and the triatomic molecule 
BCD at the classically available energy E with respect to the free fragments. Assuming that 
the harmonic approximation is valid for the intra-molecular motion of BCD, its vibrational 
energy Ey at the end of the collision reads (see appendix A) 

E v = u^ + -) + u 2 (x 2 + -) (4) 

where uj\ and uj 2 are the energy spacings between neighboring states for the two vibrational 
streching normal modes of BCD and X\ and x 2 are their related actions (since B, C and D 
are aligned, the usual bending vibration is ignored) . 

The relative translational energy Ex between A and BCD satisfies the identity 

E T = E- E v . (5) 

We call p(xi,x 2 ) the classical distribution of the actions X\ and x 2 , supposed to be 
normalized to unity. 

Additional paragraph 1: 

We shall consider the formal expressions of both the translational energy distribution 
of the final products and the one of their quantum states. However, we shall only 
represent the former distribution in the figures. We might have done the contrary, but the 
translational energy distribution is by far the most widely measured in molecular beam 
experiments. We thus believe that discussing the different ways this distribution can be 
represented in QCT studies is an important issue. In addition to that, the translational 
and internal energies being mathematically related (see Eq. ([5])), the two distributions can, 
in principle, be deduced from each other. In this section and the next one, for instance, it 
will turn out that the product state distribution is readily obtained from visual inspection 
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of the translational energy distribution. 

End of the additional paragraph 1. 



B. Purely classical translational energy distribution 

The translational energy distribution obtained from a strict application of classical 
mechanics reads 

P C (E T ) = / dxi dx 2 p(x 1 ,x 2 ) 8{Er - E + u 1 (x 1 + -) + u 2 (x 2 + -)) (6) 

(see appendix B for its derivation). Since this density has no quantum attribute, it is 
usually in bad agreement with quantum scattering and/or highly resolved experimental 
distributions, unless E is much larger than the average quantum level spacing. 

Additional paragraph 2: 

Nevertheless, this distribution has been so widely used in QCT studies that for the 
sake of completeness, we shall be considering it in this work. 

End of the additional paragraph 2. 



C. SB distributions 

A variant of the previous distribution, incorporating to some extent the idea of vibra- 
tional quantization, is as follows: 

f 1 _ 1 

Psb{E t ) = dxi dx 2 p(xi,x 2 ) 5{E T -E + ui{xi + -) + u 2 (x 2 + -)). (7) 
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We note that the only difference with respect to Eq. (^ is that the Xj's in the delta function 
have been replaced by the Xi's. 

As in the SB procedure, the bins are one unit wide, the domain of integration in Eq. ((7j) 
consists of the domains corresponding to each pair of the integer quantum numbers n% and 
n 2 . Decomposing the integral in ((7j) into a sum of integrals over these unit-sized domains, 
one gets: 



P SB (E T ) = y^ dxi dx 2 p(x 1 ,x 2 ) 5(E T - E + wi(xi + - 

ni,n 2 " , - D «i"2 



2 )+uj 2 (x2 + \)) (8) 



where D nin2 is the unit square in the plane (xi,x 2 ) centered on (711,712). The Xj's being 
equal to the rij's in D niri2 , we then arrive at 

Psb(Et) = J2 P SB(nu n 2 ) 5(E T - E + u^rn + -) + u 2 (n 2 + -)) (9) 



n 1 ,n 2 

where 



PsB{ni,n 2 )= dx 1 dx 2 p(xi,x 2 ) (10) 

" i'Tli Tin 



is recognized to be the SB population of the quantum state (%, n 2 ) 



D. GB distribution 



The GB distribution is readily found from Eqs. ([9]) and (fTPl) to be given by 

P GB (E T )=Y, P G B(n 1 ,n 2 )5(E T -E + cu 1 (n 1 + -)+u 2 {n 2 + -)) (11) 



ni,n 2 



and 



with 



PGB{n\,n 2 ) = / dx 1 dx 2 G(x 1 ,x 2 ) p(x 1 ,x 2 ) (12) 

" J-^ra no 



G(xi,x 2 ) = G(x\ - ni)G(x 2 - n 2 ). (13) 
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The even unit weight in the integrand of Eq. (TTUj) has thus been replaced by the Gaussian 
weight G(xi,X2). 

When making e tend to zero, we arrive at a distribution which we shall call "exact" in 
the following. It is of course not exact in the true quantum mechanical sense, but it is the 
best distribution we can arrive at by simple inclusion of Bohr quantization in QCTM. In 
this limit, G(u) tends to the delta-function 5(u) and Eq. (fT2|) reads 



Pe(wi,n2) = / dx\ dx2 S(xi 



-ni) S(x 2 -n 2 ) p(x 1 ,x 2 ), 



(14) 



giving immediately 
Then, Eq. (ITTj) reads 



P E (n 1 ,n 2 ) = p{n 1 ,n 2 ). 



(15) 



P E (Et) = J2 P^ H ^ 6 ( Et ~ E + Ul ( ni + o) + W 2("2 + „))• 



(16) 



rai,n 2 



Since the 1GB distribution to be derived in the next section is supposed to be an alternative 
to the GB one, the former will be systematically tested against the latter and its "exact" 
limit in the followings. 

E. 1GB distributions 

1. Statistical case 



Indirect chemical reactions involving 
intense research during the last few years 



ong-lived complexes have been the subject of 



45 



50|. 



Additional paragraph 3: 



In Phase Space Theory, the simplest statistical approach (see references 46| and 50 1 
and references therein), the final product states consistent with total energy and total 
angular momentum are equally probable. 



8 
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End of the additional paragraph 3. 

For the present system, the analogous situation corresponds to a uniform density 
p(x\,x 2 ) in the energetically available action space defined by 

E>u l {x l + -)+u 2 {x 2 + -) (17) 

and both x\ and x 2 greater than minus 1/2. This triangular domain is represented in Fig. [1] 
for E, (jj\ and u 2 kept at 5.2, 1 and 2.3 respectively ; these values have been chosen in such a 
way that six quantum states, indicated by red dots, are available (different values might have 
been chosen as well). Three forbidden quantum states are represented by dark blue dots 
and three unit squares centered on quantum states are emphasized. The two salmon ones 
correspond to available quantum states while the light blue one corresponds to a forbidden 
state. 

Throughout the present part, p(xi,x 2 ) will be simply denoted p. Its value is 

(the inverse of the area of the green triangle) inside the triangle and zero outside. 
Pc{Et) can be determined analytically by using the identity 

S(ax) = 7-rS(x), (19) 

CI 



leading to 



P c{ E T ) = l(l-^). (20) 



Additional paragraph 4: 

In principle, the translational energy distribution measured in a perfect experiment 
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would consist in a set of Dirac delta functions for the energies complementary with those of 
the allowed product quantum states. However, no experiment is "perfect". There is always 
an uncertainty in both the total energy available to the products and the measure of the 
translational energy. Consequently, the peaks are necessarily broaden. In order to take into 
account this uncertainty in the theory, we shall replace the Dirac peaks present in Eqs. (J5J), 
(ITT]) and ( 1T6|) by G(u) (see Eq. ([2])) with e = 0.05, meaning that the uncertainty on Et 
is ~ 0.1. This arbitrary value of e makes the peaks neither too narrow, nor too broad as 
compared to the total energy E. The fact that it was kept at the same value as in the GB 
procedure should not confuse the reader. Its choice was a matter of convenience, nothing 
else. 

It is worth emphasizing that there are thus two separate issues in this work: the first, 
and central one, is the use of Gaussians to deal with Bohr quantization in QCT calculations 
; the second, and minor one, is the use of Gaussians to take into account the uncertainty in 
the measurement of the translational energy. 

End of the additional paragraph 4. 

Psb{Et) and Pgb{Et) were determined by Monte-Carlo integration over x\ and x 2 , 
using Nt = 200000 points randomly chosen in the triangular domain. The corresponding 
expressions are 

P X (E T )=J2 Px(n l ,n 2 )G(E T -E + u 1 (n l + ^)+U2(n2 + \)) (21) 

n 1 ,n 2 

where X stands for SB or GB, with 

P S B(n 1 ,n 2 ) = ^^, (22) 



N' 



T 



N Dn n being the number of trajectories ending in D nin2 , and 



N D ni n 2 

P G B(n 1 ,n 2 ) = — J2 G(*i>4), (23) 



1 fc=i 
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x\ and x\ being the final actions for the k th trajectory ending in D nin2 . 

The four distributions are represented in Fig. [2] for 500 values of E t regularly distributed 
between and E. 

Pc{Et) appears to be in complete disagreement with the "benchmark" distribution 
Pe{Et), for by construction, no structure can be reproduced. Psb(Et) takes the struc- 
tures into account, but the heights of the peaks are inaccurate for the small values of Et- 
Conversely, Pqb{Et) is in excellent agreement with Pe{Et), as expected. Note that for 
these two distributions, the peaks have the same height, meaning that the populations of 
the available quantum states are all equal, in agreement with the statistical hypothesis. 

The heights of the peaks of Psb{Et) increase with E t . This can be easily understood 
from Fig. [TJ The values of the translational energies corresponding to the top of the peaks 
(see Fig. [2]) are given by 



l-<n\n2 •*-' 



( Wl (niH--)+u; 2 (n 2 H--)), (24) 



the integers n\ an n 2 being such that (ni,?^) is an allowed quantum state. These states are 
(1, 1), (3,0), (0, 1), (2,0), (1,0) and (0,0), by order of increasing translational energy. The 
straight lines defined by 

ui{xi + -) + W 2 {X2 + -) = E-E nirl2 (25) 

are represented in Fig. [TJ Clearly, the spacings between these lines exactly follow the spacings 
between the nearest i£ nin2 's (see Fig. [2]). 

Now, the heights are proportional to the SB populations, themselves proportional to the 
available areas of the unit squares centered on the quantum states. From Fig. [TJ it is clear 
that for (0, 0), this area, represented in salmon, is 1, but for (1, 1), it is only equal to ~ 0.6. 
One also guesses that for (3, 0) and (0, 1), the areas are equal to ~ 0.7 and ~ 0.9 respectively 
while for the remaining states (1,0) and (2,0), they are equal to ~ 1. This explains why 
the heights of the first three peaks of Psb{Et) (see Fig. [2]) are only ~ 60, ~ 70 and ~ 90 
percent of the height of the remaining peaks. This cannot happen with the GB distribution, 



11 



E 1GB distributions II THEORETICAL ANALYSIS OF THE 1GB PROCEDURE 

for with sufficiently thin Gaussians, their value is the same for all the quantum states. 

Having defined the system of interest and applied the different methods currently utilized 
in QCT calculations today, we are now in a position to introduce the fGB procedure. 
-PeC^i)^)) given by Eq. (TH| with p(xi,x 2 ) = p, can be rewritten as 

PE(ni,n 2 ) = p / dy 1 dy 2 <%i) 8(y 2 ) (26) 

where y\ and y 2 are two new coordinates related to x\ and x 2 by the rotation 

y\ = cosO (xi — rii) — sinO (x 2 — n 2 ) (27) 

and 

y 2 = sinO (x 1 — rii) + cosO (x 2 — n 2 ) (28) 

with 

coaB = RT^IFJ (29) 



and 



Sin0 = t 2 i 1 2M/2 - ( 30 ) 

(cof + UJ2Y' 1 



The axis y^ runs through (ra l5 n 2 ) and is parallel to the hypotenuse of the green triangle. y\ 
is thus one of the red axes represented in Fig. [lj The axis y 2 also runs through (rii, n 2 ) and 
is orthogonal to y\. These axes are represented in Fig. [3] as well as the unit square D nin2 . 
It is clear that integration over y\ and y 2 in Eq. (1261) immediately leads to Eq. (1151) . 
Let us now define the population 

P' E (ni,n 2 ) = a p dy x dy 2 S(y 2 ) (31) 

where 

a = max(cos9, sinO). (32) 

As compared with Eq. (1261) . Eq. (13 ip involves only one Dirac distribution. After a trivial 
integration with respect to y 2 , P' E (rii,n 2 ) reads 



12 



E 1GB distributions II THEORETICAL ANALYSIS OF THE 1GB PROCEDURE 



P' E (n 1 ,n 2 ) = a p / dy x . (33) 

However, one deduces from Fig. [3J corresponding to the case where cosO is greater than 
sinO {9 lower than 7r/4), that 

cos6 = - r — . (34) 

J Dnin2 *Vi 

In the same way, it can be easily shown that when cosO is smaller than sin8, 

sinO = -j. — . (35) 

Jd d y± 

a, given by Eq. ( l32l) . can thus be rewritten as 

a = — , (36) 

Jd d y± 

and Eq. ( 1331) leads to 

P' E ( ni ,n 2 ) = p. (37) 



We thus arrive at the conclusion that the P' E (ni,n2) , s are equal to the -PE(ni,ri2)'s (see Eq. 
(fT5"j) ). Consequently, P E (n 1 ,n 2 ) can be rewritten from Eqs. (I28l - (I3~T1) as 

pi \ f a a xf Ul (x 1 -n 1 )+u 2 (x 2 -n 2 )\ , . 

P E {ni,n 2 ) = a I dy 1 dy 2 d[ 7 - r - — ^m p. (38) 

JD nin2 v l w i +u 2 ) I / 

Moreover, from Eqs. (129|) . (130]) and ( 132]) and using the fact that a5(q) = 8(q/a) (deduced 
from Eq. (fl~9l)). we obtain 

PE{ni,n 2 ) = / dyi dy 2 6( p. (39) 

JD ni „ 2 v max{ujx,u 2 ) J 

Applying the GB procedure, i.e. replacing 5 by G, and going back to the Xj's, PE(ni,n 2 ) 
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finally reads 

n i \ f j j ri (u 1 (x 1 -n 1 )+uj2 {x 2 -n 2 )\ , . . s 

PiGB{ni,n 2 ) = / dxi dx 2 G[ r p(xi,x 2 ), (40) 

JD nin2 v mar(wi,w 2 ) / 

an expression formally close to the expression 

of \ f . . (u l (x l -n l )+uj 2 (x 2 -n 2 )\ 

PiGB\n 1 ,n 2 ) = dx!dx 2 G[ ■ ) p{x 1 ,x 2 ) (41) 

JD nin2 V wi+a; 2 

corresponding to the calculations of Czako and Bowman (3 
PiGB^Et) is then deduced from -Piob^i,^) by 

Pigb(Et) = J2 p iGB(n u n 2 ) 5(E T -E + u 1 (n 1 + -) + u 2 (n 2 + -)). (42) 



ni,n 2 



Like previously, the 5 function in the above expression is replaced by a Gaussian in order to 
take into account the uncertainty in the measurement of the translational energy. Moreover, 
the Monte-Carlo expression of P\GB{ n ii n 2) reads 



1712 ,, , l^. « ^ _i_, , f^fc 



P 1GB (n 1 ,n 2 ) = — ^ G( j, (43) 



1 fe=l 



with 

w = moa;(wi,W2) (44) 

according to Eq. ( I40p or 

us = us\ + uj 2 (45) 

according to Eq. ( 14"T|) . Eqs. (J2]), ( |4*3|) and (145|) correspond to Eqs. 13 and 16 in the paper 
by Czako and Bowman |37j |. The distribution obtained from Eqs. (H2"j) - (H4"|) is represented 
in Fig. H] (1GB curve), together with Pe(Et) (like previously, e was kept at 0.05 for all 
the Gaussians). The two densities are in such a good agreement that they cannot be 
distinguished. On the other hand, the distribution obtained from Eqs. (l4*2|) . (l4*3|) and (|4"5]) . 
also shown in Fig. H] (1GB' curve), has the good shape, but its norm is too large. We shall 
come back to this important normalization issue in section II. H. 

14 
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Additional paragraph 5: 

One may wonder wether the selection of y 2 to be the argument of the delta function 
in Eq. (l3Tj) is arbitrary. For instance, one might be tempted by stating in view of Fig. [3] 
that interchanging y\ and y 2 still leads to the equality between Pg (711,712) and P^ (711,712). 
This is indeed true for the less excited states (0,0), (1,0), (2,0) and (0,1), but not for 
the most excited states (3,0) and (1, 1) (see Fig. [TJ. As a matter of fact, the upper limit 
of the green triangle would limit the integration along the new y\ axis (corresponding to 
the 2/2 axis in Fig. [3]) in such a way that for (3,0) and (1,1), Pg (711,712) would be lower 
than Pg (711,712). On the other hand, slightly varying 9 about the value defined by Eqs. 
( 129]) and (130]) will preserve the equality between Pg (711,712) and P' E (ni,n 2 ) provided that 
the hypotenuse of the green triangle is not too close to the most excited states. Strictly 
speaking, the value of 9 defined by Eqs. ( 129]) and ( 130]) is thus not the only satisfying one. 
Nevertheless, it is the only one for which the equality between Pg(ni, 71-2) and P' E (ni, 712) is 
systematically satisfied, no matter how close to the most excited states the hypotenuse is. 
This makes it superior to any other one. 

End of additional paragraph 5. 



2. Non statistical case 

What about non statistical situations ? Given that p(xi,x 2 ) is non uniform in the 
energetically allowed triangle, Eq. ( 133]) reads 

P' E (n 1 ,n 2 ) = a / dy x p(x x , x 2 ) . (46) 

Recall that in the statistical case, the 1GB procedure was justified by the fact that Pg(rii, 712) 
is equal to P E (ni, n 2 ). 

In the present case, P E (ni,n2) is still a reasonable approximation of Pe(tii,712) provided 

15 
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that the variation of p{x\,x-i) along the 2/1-axis is sufficiently smooth. If so, it is indeed 
clear using Eq. ( 136|) that 

f Id„ „ d 1h P^i)^) 

a/ dy 1 p{x 1 ,x 2 ) = — ^ & p(rii,n 2 ), (47) 

the strict equality occuring when p(xi,X2) varies linearly along the ^-axis (and is non zero 
within D nin2 ). From Eq. (1T5|) . we then arrive at the conclusion that P l E {n\ 1 n2) is roughly 
equal to P^rii,-^). 

Such a smooth variation of the density in the action space was observed several times 
in the case of three-atom exchange reactions, the only difference being that the previously 
introduced x and j actions play the role of X\ and x<i- Fig. 4 in reference 20| is a clear 



illustration of this statement in the case of the direct reaction 0( 3 P)+HC1 — > OH+Cl( 2 P) 
The energetically available area in the (x, j) plane is given by 



E > oo{x + -) + Bj 2 . (48) 



Its upper limit is thus found to be a curved, instead of straight, line. Moreover, the lines 
equivalent to the six straight lines in Fig. [T] are also curved. Along these lines (not drawn 
in Fig. 4 of reference [20]), the density p(x,j) evolves rather smoothly. In other words, 
despite the fact that the distribution of the translational energy is not statistical, the intra- 
molecular vibrational redistribution (IVR) in the strong coupling region tends to distribute 
in a relatively democratic way the rest of the energy among the vibrational and rotational 
degrees-of- freedom. 

One expects a similar redistribution will also take place among the vibrational modes, 
due to their couplings. As shown in section IV, this is at least the case for OH+D2 which, 
as 0( 3 P)+HC1, is a direct process. 

F. Collisional system involving three vibrational modes 

We now consider the collinear inelastic collision between atom A and the tetra-atomic 
molecule BCDE involving three vibrational normal modes. Following a reasoning analogous 
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to the one developed in the previous subsection, we arrived after some steps of algebra to 
an expression exact (with a Gaussian of zero width) in the statistical limit (Pe{Et) being 
still considered as the "exact" density). With n = (m, 712,^3), x = (x\,X2, X3), D n the unit 
cube centered on n, and uj^, uji and uo m deduced from any cyclic permutation of uji, 002 and 
W3, this expression is 

3 

Pigb{E t ) = J2 p iGB(n) S\Er-E + J>^ + -)) (49) 

n «=1 

with 

P 1GB (n) = f dx G ( P=^ {Xt ~ ni) ) p(x), (50) 



uj being defined as 



^ = 7 w (51) 



if Uk < oji + u m for the three cyclic permutations, or 

uj = maxiujk, u>i, cu m ) (52) 

if ^fc > Ui + bj m for only one of the permutations. In the second case, the formulation (see 
Eqs. (150|) and ( 152]) ) is a straightforward extension of Eq. ( |4"0]) . 

We shall retain from the above developments that for three vibrational modes, the for- 
mulation of uj is not unique. It depends on the values of the vibrational frequencies Ui, U2 
and CU3, contrary to the formulation for two vibrational modes. 

Hence, if one does not use the correct expression of u, one does not find the correct 
populations. However, the wrong populations turn out to be proportional to the correct 
ones. The proof is straightforward: if we call u c the correct value of uj and uj w the wrong one, 
the wrong populations are found from Eqs. ( 1501 (with G replaced by the Dirac distribution) 
and ( TT9l) to be equal to uj w /uj c times the correct populations. This explains why the peaks 
of the 1GB' distribution in Fig. H] (given by Eqs. (T4TT) and (142]) ) are (ui + UJ2)/max(uJi, UJ2) 
higher than the peaks of the 1GB distribution (given by Eqs. ([4*0]) and ([4*2]) ). 

We did not extend the above developments in the case of systems involving more than 
three vibrational modes for the mathematical developments became very tedious. Therefore, 
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we do not know the analytical expression of u making the 1GB distribution in close agree- 
ment with the "exact" or GB distribution in the statistical limit. However, we go round this 
difficulty in section II. H. 

G. General collisions 

The extension of Eq. (jSJ) to a three-dimensional collision involving iV vibrational modes 



is 



P C (E T ) = J dE R P C (E T , E R ) (53) 



where 

N 

P C (E T ,E R ) = / dxp(x,E fl ) 5(Er -{E- Er) +Y,u i {x i + -)), (54) 

x = (xi, ...,Xn), E r is the final product rotational energy and p(x,Er) is the distribution 
of x and Er. 

For a given value of Er, Pc{Et,Er) appears to be formally identical to Pc{Et) in 
Eq. ([6]) and consequently, all the developments following Eq. ([6]) could be repeated here 
identically. The main conclusion of this section is thus the same as before, i.e., the 1GB 
procedure leads to nearly the same conclusions as the GB procedure provided that the 
variation of p(x, Er) in any plane parallel to the plane 

N 1 

E = E R + J2 UJ n(x i + -) (55) 

i=i 

is sufficiently smooth. As stated before, however, we do not know the expressions analogous 
to Eq. gOD and Eqs. (15t)|) - (lH2]) for iV larger than 3. 

H. Normalization procedure in realistic calculations 

The Monte-Carlo expression of the 1GB populations of the final product quantum states 
n = (m, ■■■,n N ) is given by 

P 1CB(n) = J- EG ( ^^.-) ) (56) 

1 fc=l 
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where N? is the total number of trajectories run and N n is the number of trajectories ending 
in the product channel with x = (xi, ...,Xn) pointing in D n , the iV-dimensional unit cube 
centered on n. 

As seen before, however, we do not know in the general case the expression of a; leading to 
1GB distributions in close agreement with the GB ones in the statistical limit. Nevertheless, 
the former are proportional to the latter. 

One might thus think about re- normalizing 1GB distributions so as to give them the 
GB norms. But GB norms have no reason to be exactly equal to one, so it is preferable to 
directly re-normalize 1GB distributions to unity. The corresponding expression is 

PigbH = * =1 )_„ " .{ (57) 

where u may be kept at the maximum of the frequencies or their sum. The final result is 
not expected to depend significantly on this choice, for it will affect both the numerator and 
the denominator of Eq. (157)) in nearly the same way. This point is illustrated in section IV 
in the case of the reaction OH + D 2 . Note that in the denominator, the sum is over the 
whole set of computed trajectories, be they reactive or not. The various quantities in the 
argument of the Gaussian are thus either those of the products and those of the reformed 
reagents. 

Special care should however be taken with processes involving a large amount of 
vibrationally elastic non reactive trajectories. Ion-molecule reactions are a typical example. 
For such processes, an alternative to Eq. (157)) is 



mm 
iWn) 



N n ,E^G( 



N n WEJIi^O^-rcQ 



E m min [jV m ,E£i G( 



J2^Li^t( x i- m i) 



(58) 



where the sum over m in the denominator involves the whole set of final vibrational states, 
i.e., those of the products as well as those of the reformed reactants. Eq. (158)) is a compact 
form of Eqs. (15) and (16) of reference 23. Eq. (4) of reference 34 may be a second 
alternative. 
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III. NON STATISTICAL TEST CASE FOR THREE VIBRATIONAL MODES 

We still consider the collinear inelastic collision between atom A and the tetra-atomic 
molecule BCDE involving three harmonic normal modes. E, Ui, u 2 and U3 are respectively 
kept at 15, 1, 1.7 and 2.9. We also consider the non statistical Gaussian density p(x\, x 2 , £3) 
given by 

p(x 1 ,x 2 ,x 3 ) =T$ =1 G(xi-Xi) (59) 

with x\ = 2.2, x 2 = 1.3 and 2° = 0.7, e being kept at 0.8 in G(x 1 — xf), 1.4 in G(x 2 — x 2 ) 
and 0.4 in G(x^ — x®). 

The resulting distributions P C (E T ), P SB (E T ), P GB (E T ), P 1GB (E T ) and Pe(E t ), given by 
expressions similar to those of the previous section with one more dimension, are represented 
in Fig. [5j The details of the calculations are exactly the same as in section II.D.l, the only 
difference being that all the distributions were numerically re-normalized to unity. 

Like in the previous statistical example, Pc(Et) is in poor agreement with Pe{Et)- On 
the other hand, Pgb(Et) is in very good agreement with Pe{Et) ; given the large number 
of points considered in the Monte-Carlo integration, this is an expected result despite the 
already "large" number of vibrational modes involved in the collision. 

Interestingly, Piqb{Et) is even in slightly better accord with Pe(Ej-) than Pgb{Et) 
when looking at the details. Such a high level of agreement despite the non-statistical 
nature of the present process is pleasing. It supports the statement of subsection II. E. 2 that 
for a sufficiently smooth distribution of the final actions, the 1GB procedure represents an 
accurate alternative to the GB one. 

Last but not least, Pgb{Et) does also a good job, though the heights of the peaks 
corresponding to the largest energies tend to be overestimated. 

IV. THE SIMPLEST POLYATOMIC REACTION OH+D 2 — ► HOD+D 

This process has been the subject of intense research, both experimentally and theoreti- 



cally 



38 



44 



5l|. Its mechanism has been well established as being direct, with the products 



preferentially backward scattered, suggestive of a rebound mechanism. The product trans- 
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lational energy distribution was experimentally measured by Alagia et al. 38j, and by Davis 



and co-workers four years later J40|. While Alagia et al. found a broad and single-peaked 
distribution, Davis et al. found a better resolved distribution involving three peaks corre- 
sponding to the HOD vibrational states (ni,n 2 ,n 3 ) = (0,1,0), (0,2,0) and (0,1,1). n 1; n 2 
and ri3 are the OH stretching, bending and OD stretching quantum numbers, respectively. 
This distribution is represented in the top panel of Fig. [6l 

In reference 41, 1 000 000 trajectories were run on the Ochoa-Clary (OC) potential 
energy surface (PES) |51| using the VENUS96 code. Initial conditions were selected to 
reproduce the experiment of Davis and co-workers, with a collision energy of 6.6 kcal mol^ 1 
and the reactants in their vibrational ground states. The number of reactive trajectories 
was found equal to Np = 10837. At the end of each reactive trajectory, the vibrational 
actions of the triatomic HOD product were calculated using the recent normal mode analysis 
algorithm 



52| . The latter includes anharmonicity and Coriolis-coupling terms, and yields 
results similar to those obtained by means of the widely used fast Fourier transform approach 
[53|, but at a lower computational cost. 

The different distributions previously considered are calculated as follows. Formally, the 
purely classical translational energy distribution is given by Eqs. (153]) and (154"|) with N 
equal 3. Stricto-sensu, its Monte-Carlo expression is 

P C {E T ) = W Y, S &r ~ E T ) (60) 

T fe=i 

where E T is the final translational energy for the k th trajectory. This energy satisfies the 
relation 

3 1 

E k T = E - Yl <*(*? + o) - E ^ ( 61 ) 

the x\ 's and E\ being the final vibrational actions and rotational energy for the k th trajec- 
tory. 

Here, we shall not replace the Dirac distribution in the previous sum by a Gaussian and 
calculate it for fixed values of Et- Instead, we divide the available range of energy [0, E] in 



N r boxes [(i — l)E/N r ,iE/N r ], i = l,N r , and integrate Pc(Et) over the boxes. This leads 
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to the N r populations 



p - = W < 62 > 



i = l,N r , where Ni is the number of trajectories for which the final translational energy 
belongs to the i th box. 

We note that applying Eq. ( J60l) does only require the calculation of Ej>, not of x\, x%, 
X3 and E\. On the other hand, the four last quantities are necessary for the calculation 
of P$b(Et). This distribution is indeed calculated in the same way as Pc{Et), the only 
difference being that the translational energy for the k th trajectory is now given by 



Ulrp 



3 1 

E - J2 <*(*? + o) - &R (63) 



8=1 



instead of Eq. ( ISTj) . This difference is similar to the one between Eqs. (Q and ([7|). 
For P GB (E T ), N { is replaced in Eq. ([62]) by 



Ni 



A?* = £llLi <?(*}-*?), (64) 



ft=l 



the sum being performed over the trajectories for which the final translational energy ac- 
cording to Eq. ( ISTT) belongs to the i th box. 
For Picb{Et), Ni is replaced by 



N ? B = E G ( Et=1 ^ ( ^~^ ) , (65) 



fc=i 

the sum being performed over the same trajectories as above. 

Finally, Psb{Et), Pgb{Et) and Pigb(Et) were re-normalized to one. e was kept at 
0.05 for Pqb(Et), and 0.01 for P\.gb(Et)- w was identified with the largest frequency. 
The distributions are displayed in Fig. [6j We also kept 00 at the sum of the frequencies, 
following Czako and Bowman 



37J, but due to the re- normalization, this left the distribution 



unchanged. 
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Contrary to the purely classical distribution which has a bell shape and does not re 



produce the vibrational structures observed experimentally [43l. I44J . the SB, GB and 1GB 
distributions reproduce quite satisfyingly the two structures due to the (0,1,0) and (0,2,0) vi- 
brational states. On the other hand, the third structure, due to the (0,1,1) state, is strongly 
underestimated by all the treatments. Comparison with exact quantum scattering calcu- 
lations, certainly possible in a near future, will tell if the previous disagreement is due to 
possible inaccuracies of the OC-PES or to the present classical descriptions. 

The SB procedure does a good job, though it overestimates the contribution of the small 
translational energies to the (0,2,0) peak. 

The GB distribution of Fig. [6] involves strong fluctuations, contrary to the same density 
represented in Fig. 2 of reference 41. The reason is that in the present work, we did not 



use the smoothing procedure previously considered 44j (two Gaussian functions were used 
to fit the left and right-hand side of each vibrational contribution). We did it on purpose, 
in order to illustrate the fact that with ~ 11000 reactive trajectories and three vibrational 
modes, the usual GB procedure generates quite noisy curves. On the other hand, the 1GB 
distribution is much better converged and one guesses that it represents the curve one would 
obtain from smoothing the GB curve. 

Owing to the fact that the OH+D2 reaction is a direct process, the present results are 
quite encouraging for future applications of the 1GB procedure to polyatomic reactions, do 
they involve a long-lived complex or not. 

V. CONCLUSION 

In the recent years, many processes have been studied by the quasi-classical trajectory 
method (QCTM) within the Gaussian binning (GB) procedure. In most studies, the 
population of the final product quantum state n = (n\, ...,Un), N being the number of 
quantized degrees of freedom (DOF), was approximated by 

^ = j- E nf =1 G(x1 - m) (66) 

T k=i 

instead of the usual expression 

p » = w T (67) 
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used in the standard binning (SB) procedure (or histogram method). N? is the total number 
of trajectories run, x = (a?i, ..., xn) is the final action state, N n is the number of trajectories 
ending in the product channel with x pointing in the N- dimensional unit cube centered on 
n and G is a Gaussian normalized to unity, with a full width at half maximum of ~ 10 
percent. 

Since most processes studied so far by GB-QCTM were triatomic reactions, one single 
vibrational DOF was quantized, meaning that the Gaussian product in Eq. f)66|) reduced to 
one term only. As about 10 percent of the total amount of reactive trajectories do actually 
contribute to the product populations, 10 times more trajectories had to be run for keeping 
with the same level of convergence of the predictions as compared with SB-QCTM. 

Nowadays, however, more and more processes under scrutiny involve more than one 
vibrational mode. For instance, the reaction OH+D 2 — > HOD+D involves three modes 
while for the reaction F+CH4 — y FH+CH3, this number is seven. Consequently, GB- 
QCTM requires one thousand more trajectories than SB-QCTM for the first process and 
ten millions more for the second ! It is thus quite clear that as such, Eq. (166]) has no future 
in the area of polyatomic reaction dynamics. 

This is why Czako and Bowman [37] recently proposed to "quantize" the total vibrational 
energy instead of the vibrational actions, introducing the expression 






■4£ ^?^ ) <«» 



1 fc=l 



where Wj is the i th normal mode frequency and 

N 

J>. (69) 



N 



The key feature of this ad-hoc quantization as compared to the previous one-Gaussian- 
for-one-mode approach is that only one Gaussian function is used whatever the number of 
vibrational DOF of the system, a huge amount of computational time being therefore saved. 
We called it the 1GB procedure. 

The conclusions of the present paper are as follows: 
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1) For a statistical collision involving two product vibrational modes, the 1GB pro- 
cedure is strictly equivalent to the GB procedure provided that u is identified with the 
maximum of the Wj's instead of their sum. 

2) For a statistical collision involving three product vibrational modes, the 1GB pro- 
cedure is strictly equivalent to the GB procedure provided that u is kept at the maximum 
frequency in part of the frequency space, and a more complex expression (see Eq. (1511) ) in 
the remaining part. 

3) For the previous processes and a non statistical but sufficiently smooth distribu- 
tion in the action space, the 1GB procedure leads to results in satisfying agreement with 
the GB ones. 

4) Finding the expression of u for any realistic process involving more than three 
product vibrational modes requires heavy mathematical developments that we did not 
perform. However, one may go round this difficulty by re-normalizing the product state 
populations. In such a case, u can be indifferently kept at the maximum of the Wj's or 
their sum. Special care should however be taken with processes involving a large amount 
of vibrationally elastic non reactive trajectories, like ion-molecule reactions. The meth- 
ods proposed in reference 23 (leading to Eq. ( |58|) of the present work) or 34 can then be used. 

5) The 1GB procedure leads to results in good agreement with the GB one for (a) a 
non statistical test case involving three vibrational modes and (b) the prototype four-atom 
reaction OH+D 2 — > HOD+D. 

In conclusion, the 1GB procedure might be of great interest for future classical sim- 
ulations of polyatomic chemical reaction dynamics in the highly quantum mechanical 
situation where only a few product vibrational states are available. 
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Appendix 

Consider (i) the iV-dimensional space T = (xi, ...,Xn), (ii) a given distribution p(T), 
normalized to unity, of the position in the previous space and (iii) the quantity Q' depending 
on r according to 

Q> = /(T). (B.l) 

The probability that Q' is lower than a given value Q is given by 

II(Q) = j dT p(T)Q(Q -Q') = J dT p(T)Q(Q - f(T)) (B.2) 

where Q(x) is the Heaviside function, equal to for x negative and 1 in the contrary case. 
Q(Q — Q') ensures that integration with respect to T is made over the volume such that Q 
minus Q' is positive, i.e., Q' is lower than Q. 

If P(Q) is the density of probability that Q' takes the value Q, P(Q)dQ is the probability 
that Q' belongs to the range [Q, Q + dQ\. We then have 

P{Q)dQ = U(Q + dQ) - n(Q) (B.3) 

that is, 

P(Q) ^ %. (B.4) 



From Eq. (1B.2J) . we finally arrive at 

P(Q) = J dTp(T)6(Q-f(T)), (B.5) 

as the Dirac distribution 5(x) is the first derivative of Q(x). Eq. (jBJ) is straightforwardly 
obtained from Eqs. ( IB. 51) . (jlj) and (jSJ). 
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Figures captions 

Fig. [JJ The action space denned by Eq. ( |17p and both X\ and x 2 greater than 
minus 1/2 is represented here by a green triangle while the available quantum states are 
represented by six red dots. The distribution of the action pair {x\,x-i) is uniform in the 
triangle. Each dot lies along a straight-line corresponding to a given translational energy 
(see Eq. (125]) ). Three forbidden quantum states are represented by dark blue dots and 
three unit squares centered on quantum states are emphasized. The two salmon ones cor- 
respond to available quantum states while the light blue one corresponds to a forbidden state. 

Fig. |2) Translational energy distributions corresponding to the statistical distribution in 
the green triangle of Fig. [JJ The curves are labeled by the subscript of their mathematical 
symbols (see text). The distributions are not normalized to unity. 

Fig. [3) Drawing shedding light on the derivation of Eq. ([4*01 . cosO clearly appears to be 
the reciprocal of the section of the new coordinate axis y\ being within the red square (see 

Eq. AMD). 

Fig. HJ Translational energy distributions corresponding to the statistical distribution in 
the green triangle of Fig. [TJ The curves are labeled by the subscript of their mathematical 
symbols (see text). The distributions are not normalized to unity. 

Fig. |5j Translational energy distributions corresponding to the non statistical distri- 
bution given by Eq. (159|) . The curves are labeled by the subscript of their mathematical 
symbols (see text). The distributions are not normalized to unity. 

Fig. El Translational energy distributions in the products of the reaction OH+D2 
— > HOD+D studied at the conditions of the group of Davis 



40] . The top curve is the 



experimental distribution while the remaining curves are labeled by the subscript of their 
mathematical symbols (see text). The distributions are normalized to unity. 
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